Read Pressure gradient wave and perform Fourier analysis

SetDirectory["/Users/nikosstergiopulos/DeskTop"]

filename = "Pressure Gradient CC.txt"

database = ReadList[filename,{Number,Number,Number}];

t=database[[Range[Length[database]],1]];
pz=database[[Range[Length[database]],2]];
Qexp=database[[Range[Length[database]],3]];

np=Length[pz]

ListPlot[-pz,Joined->True,PlotRange->All]

period=t[[np]];

freq=1/period

pzf = N[InverseFourier[pz]/Sqrt[np]];

ListPlot[Abs[pzf],PlotRange->{0,200}]

(* Select number of harmonics nhar *)

nhar=21;

pzfs=Flatten[{pzf[[1]],Table[2 pzf[[i]],{i,2,nhar}]}];

Rheological parameters & geometry

rho=1060;
r0=0.0055;
mu=0.004;

omega[i_]:=2*Pi*freq*i
alpha[i_]:=Sqrt[omega[i]*rho*r0^2/mu]

omega[1]

alpha[1]

Calculate the velocity profile using Womersley's theory

Get the velocity profile for each harmonic vharm

vharm[r_,t_,i_]:=Re[pzfs[[i+1]]/(I*rho*omega[i])*(1-BesselJ[0,I^(3/2)*alpha[i]*r/r0]/BesselJ[0,I^(3/2)*alpha[i]])*Exp[I*omega[i]*t]]

Synthesize velocity from mean flow (Poiseuille) and all harmonics (Womersley)

v0[r_]:=pzfs[[1]]/(4*mu)*(r0^2-r^2)

v0[0]

v[r_,t_]:=v0[r]+(∑_(i=1)^(nhar-1))vharm[r,t,i]

Plot[v[r,period/4],{r,-r0,r0},PlotRange->{-0.2,0.4}]

Animate[Plot[v[r,t],{r,-r0,r0},PlotRange->{-0.2,0.4}],{t,0,period}]

Calculate flow rate

Q[t_]:=Re[NIntegrate[v[r,t]*2*Pi*r,{r,0,r0}]]

Q[period]

Plot[Q[t]*10^6,{t,0,period},PlotRange->All]

ListPlot[Qexp,Joined->True]